Read in network information for the L/R side

# Get network information for the L/R side 
path = paste0("../../Processed_FunctionalData/",params$subj_name,"/")
path_vec = c(path)

edge_time_mat_list = list()
locs_mat_list = list()
mnx_vec_list = list()
dFF_mat_list = list()
iso_inds_list = list()
non_iso_inds_list = list()


for(m in 1:length(path_vec)){ 
  path = path_vec[m]
  
  ### Read information from data
  edge_time_mat = as.matrix(read.csv(paste(path, '/EdgeTime.csv', sep='')))
  edge_time_mat = edge_time_mat[,-1]
  
  avai.inds = as.matrix(read.csv(paste(path,'/AvaiNeurons.csv',sep='')))
  avai.inds = avai.inds[,-1];
  
  locs.all = as.matrix(read.csv(paste(path,'/locs_all.csv',sep='')))
  locs.all = locs.all[,-1]
  locs_mat = locs.all[avai.inds,]
  
  mnx.all = as.matrix(read.csv(paste(path,'/mnx.csv',sep='')))
  mnx.all = mnx.all[,-1]
  mnx_vec = mnx.all[avai.inds]
  
  dFF.all = as.matrix(fread(paste(path,'/dFF.csv',sep='')))
  dFF.all = dFF.all[,-1]
  dFF_mat = dFF.all[avai.inds,]
  
  ### Split each zebrafish into left and right
  inds_L = which(locs_mat[,2]<0)
  inds_R = which(locs_mat[,2]>0)
  
  edge_time_mat_L = edge_time_mat[inds_L, inds_L]
  edge_time_mat_R = edge_time_mat[inds_R, inds_R]

  locs_mat_L = locs_mat[inds_L,]
  locs_mat_R = locs_mat[inds_R,]
  
  mnx_vec_L = mnx_vec[inds_L]
  mnx_vec_R = mnx_vec[inds_R]
  
  dFF_mat_L = dFF_mat[inds_L,]
  dFF_mat_R = dFF_mat[inds_R,]
  
  
  iso_inds_L = which(rowSums(edge_time_mat_L<Inf)<=0)
  iso_inds_R = which(rowSums(edge_time_mat_R<Inf)<=0)
  iso_inds = which(rowSums(edge_time_mat<Inf)<=0)
  
  non_iso_inds_L = setdiff(1:nrow(edge_time_mat_L), iso_inds_L)
  non_iso_inds_R = setdiff(1:nrow(edge_time_mat_R), iso_inds_R)
  non_iso_inds = setdiff(1:nrow(edge_time_mat), iso_inds)
  
  edge_time_mat_L = edge_time_mat_L[non_iso_inds_L, non_iso_inds_L]
  edge_time_mat_R = edge_time_mat_R[non_iso_inds_R, non_iso_inds_R]
  
  edge_time_mat_list[m] = list(edge_time_mat)
  locs_mat_list[m] = list(locs_mat)
  mnx_vec_list[m] = list(mnx_vec)
  dFF_mat_list[m] = list(dFF_mat)
  
  iso_inds_list[m] = list(iso_inds)
  non_iso_inds_list[m] = list(non_iso_inds)
  
}

Read in data analysis results

### Read in data analysis results
data_folder = paste0(params$data_folder, params$subj_name, "/") 
path_vec = list.files(data_folder, full.names = TRUE, recursive = TRUE)


res_list = list()

for(m in 1:1){ 
  path = path_vec[m]
  ### Read in results from data
  load(path)
  res_list[[m]] = res
}


### Collect from all subjects: clusters_list, center_pdf_array, v_vec
clusters_list = lapply(res_list, function(res) res$clusters_list)
center_pdf_array_list = lapply(res_list, function(res) res$center_pdf_array)
v_vec_list = lapply(res_list, function(res) res$v_vec)
t_vec = res$t_vec
n0_mat_list = lapply(v_vec_list, function(v_vec) n0_vec2mat(n0_vec = v_vec/(t_vec[2]-t_vec[1])))

# ### Fix spike intensity 
# for (ind in 1:length(center_pdf_array_list)) {
#   if (dim(center_pdf_array_list[[ind]])[3]<length(t_vec)) {
#     center_pdf_array_list[[ind]] = get_center_pdf_array_v2(edge_time_mat_list = edge_time_mat_list[ind], clusters_list = clusters_list[ind], n0_mat_list = n0_mat_list[ind], t_vec = t_vec)
# }
# }

Visualize estimated intensities

if (!exists("t_vec")) {
  t_vec = seq(0,336,1)
}
### Get connecting probabilities
order_clus_1 = order(rowSums(apply(center_pdf_array_list[[1]], 
      MARGIN = c(1,2),
      function(vec)sum(vec*(t_vec[2]-t_vec[1])))), 
      decreasing = TRUE)

### Visualize estimated connecting intensities 
clus_size_list = lapply(clusters_list, function(clusters)sapply(clusters,length))

g1 = plot_pdf_array_v2(pdf_array_list = center_pdf_array_list[[1]][order_clus_1,order_clus_1, ,drop=FALSE], 
                        pdf_true_array = center_pdf_array_list[[1]][order_clus_1,order_clus_1, ,drop=FALSE],
                        clus_size_vec = clus_size_list[[1]][order_clus_1],
                        t_vec = t_vec, x_lim = c(), y_lim = c(0,max(unlist(center_pdf_array_list[[1]]))))
g1

Visualize estimated time shifts

### Get memberships for ALL nodes
membership_list = lapply(clusters_list, clus2mem)
membership_list[[1]] = clus2mem(clusters_list[[1]][order_clus_1])


# for (m in 1:length(membership_list)) {
#   N_node = nrow(locs_mat_list[[m]])
#   mem_tmp = rep(0,N_node)
#   mem_tmp[non_iso_inds_list[[m]]] = membership_list[[m]]
#   membership_list[[m]] = mem_tmp
# }
### Plot time shifts distribution
# data = tibble(v_vec=unlist(v_vec_list), 
#               membership=c(membership_list[[1]][non_iso_inds_list[[1]]], membership_list[[2]][non_iso_inds_list[[2]]]), 
#               subj=rep(1:length(v_vec_list), times=sapply(v_vec_list,length)))

data = tibble(v_vec=unlist(v_vec_list), 
              membership=c(membership_list[[1]]), 
              subj=rep(1:length(v_vec_list), times=sapply(v_vec_list,length)))

g4 <- data %>%
  ggplot(aes(x=v_vec,group=as.factor(membership),
                     color=as.factor(membership)))+
  geom_density( aes(), alpha=0.6, position = 'dodge') +
  scale_fill_manual(values=palette()[1+1:3]) +
  theme_bw()+
  theme(legend.position = "none") +
  # facet_wrap(~subj) +
  scale_x_continuous() +
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank(), axis.title.y = element_blank()) +
  xlab('Time shifts (min)')

print(g4)
## Warning: Width not defined. Set with `position_dodge(width = ?)`

### Get Wan's result
dat=readMat(paste('../../Raw_data/to_shizhe/','Leader_',substr(params$subj_name, start=nchar(params$subj_name)-7, stop=nchar(params$subj_name)),'.mat',sep=''))
activeTime_vec = dat$activeTime * 60
patternTime_vec = dat$patternTime * 60

# ### Let time shift for non-active nodes be NA
# for (m in 1:length(v_vec_list)) {
#   N_node = nrow(locs_mat_list[[m]])
#   v_vec_tmp = rep(NA,N_node)
#   v_vec_tmp[non_iso_inds_list[[m]]] = v_vec_list[[m]]
#   v_vec_list[[m]] = v_vec_tmp
# }

ggplot(data = data.frame(activeTime_vec=activeTime_vec,
                         v_vec=unlist(v_vec_list),
                         mem=as_factor(unlist(membership_list))
                         )) + 
  geom_point(aes(x=activeTime_vec, y=v_vec, color=mem)) +
  # scale_color_manual(values=palette()[c(8,2:6)]) +
  scale_color_manual(values=palette()[c(2:6)]) +
  xlab('Activation time by Wan (pp. e7)') +
  ylab('Our estimation')
## Warning: Removed 10 rows containing missing values (geom_point).

ggplot(data = data.frame(activeTime_vec=activeTime_vec,
                         patternTime_vec=patternTime_vec,
                         v_vec=unlist(v_vec_list),
                         mem=as_factor(unlist(membership_list))
                         )) + 
  geom_point(aes(x=patternTime_vec, y=v_vec, color=mem)) +
  # scale_color_manual(values=palette()[c(8,2:6)]) +
  scale_color_manual(values=palette()[c(2:6)]) +
  xlab('Patterning time by Wan (pp. e7)') +
  ylab('Our estimation')
## Warning: Removed 18 rows containing missing values (geom_point).

{r} # ind_vec = which(unlist(v_vec_list)!=0 & activeTime_vec==1) # ind = ind_vec[1] # # for (i in 1:min(length(ind_vec),3)) { # ind = ind_vec[i] # data_trace = tibble(dFF=dFF_mat[ind,],time=(1:ncol(dFF_mat)*0.25)/60) # g = ggplot(data_trace, aes(x=time, y=dFF)) + # geom_line(alpha=0.3) + # # xlim(unlist(v_vec_list)[ind]+c(-1,1)) + # geom_vline(xintercept = unlist(v_vec_list)[ind]) # print(g) # } # # #

Visualize estimated clusters with spatial location

### Spatial location  ----

data = as.data.frame(reduce(locs_mat_list,rbind))
colnames(data) = c('AP','LR','DV')
data = cbind(data, 
             membership=as.factor(unlist(membership_list)),
             subj=rep(1:length(membership_list),
                      times=sapply(membership_list,length)))

# lims = data.frame(coord=c("AP","DV","LR"), xmin=c(-0.5, 0, -2),xmax=c(9,140,3))
g<- data %>% 
  pivot_longer(cols = c('AP','LR','DV'), names_to = "coord", values_to = "loc") %>%
  ggplot(aes(x=loc,
                     group=membership,
                     color=membership))+
  geom_density( aes(), alpha=0.6, position = 'dodge') +
  scale_color_manual(values=palette()[c(8,2:6)]) +
  theme_bw()+
  theme(legend.position = "none") +
  facet_wrap(~coord, scales='free',ncol=1) +
  xlab('') +
  # scale_x_continuous(limits = c(-2, 6))+
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank(), axis.title.y = element_blank()) +
  scale_y_continuous(position = "right")
print(g)
## Warning: Width not defined. Set with `position_dodge(width = ?)`

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:igraph':
## 
##     groups
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
fig = plot_ly(data, x = ~LR, y = ~AP, z = ~DV, color = ~membership, 
              # colors = c(palette()[c(8,2:6)]),
              colors = c(`0`=palette()[8], 
                         `1`=palette()[2],
                         `2`=palette()[3],
                         `3`=palette()[4],
                         `4`=palette()[5],
                         `5`=palette()[6]),
              type = "scatter3d")
fig
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Visualize estimated clusters with mnx+/- indicator

### mnx decomposition  ----

data = tibble(subj=rep(1:length(membership_list),sapply(membership_list,length)),
              membership=unlist(membership_list),
              mnx=as.factor(unlist(mnx_vec_list)), 
              islet=NA) 

data_2 = data %>% 
  mutate(mnx=factor(mnx,levels=c(1,0)), islet=as.factor(as.numeric(islet)), 
         membership=factor(membership, levels=c(0,max(membership):1))) 

data_2_left = data_2 %>% filter(subj==1) %>% mutate(membership=fct_drop(membership))
data_2_right = data_2 %>% filter(subj==2) %>% mutate(membership=fct_drop(membership))

fill_color = c(`0`=palette()[8],
               `1`=palette()[2],
               `2`=palette()[3],
               `3`=palette()[4],
               `4`=palette()[5],
               `5`=palette()[6])
### Left spine
g1 = data_2_left %>%
  ggplot(mapping = aes(fill=membership,x=mnx)) +
  scale_fill_manual(values=fill_color) +
  geom_bar(position="fill") +
  geom_text(data = data_2_left %>% 
              count(mnx, membership) %>% group_by(mnx) %>% 
              summarise(perc=round(n/sum(n),3), membership=membership) %>% 
              arrange(desc(membership), .by_group = TRUE) ,
            aes(x = mnx, label = scales::percent(perc), y=perc, fill=NULL), 
            position = position_fill(0.5), size=2.5)+
  scale_x_discrete(labels=c("mnx+","mnx-"))+
  xlab("")+
  theme_bw()+
  theme(legend.position = 'none', 
        panel.grid  = element_blank(),
        # axis.text = element_blank(),
        axis.text.y = element_blank(), 
        axis.title = element_blank(),  
        axis.title.y = element_blank(), 
        axis.ticks.y = element_blank())
## `summarise()` has grouped output by 'mnx'. You can override using the `.groups` argument.
g1